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Abstract 

We investigate strained heteroepitaxial crystal growth in the 
framework of a simplifying (l+l)-dimensional model by use of off- 
lattice Kinetic Monte Carlo simulations. Our modified Lennard- 
Jones system displays the so-called Stranski-Krastanov growth 
mode: initial pseudomorphic growth ends by the sudden appearance 
of strain induced multilayer islands upon a persisting wetting layer. 

1 Introduction 

In addition to its technological relevance, epitaxial crystal growth is highly 
attractive from a theoretical point of view. It offers many challenging open 
questions and provides a workshop in which to develop novel methods for 
the modeling and simulation of non-equilibrium systems, in general. 

In particular, strained heteroepitaxial crystal growth attracts significant 
interest as a promising technique for the production of, for instance, high 
quality semiconductor films. Recent overviews of experimental and the- 
oretical investigations can be found in, e.g., JH El El El. A particular at- 
tractive aspect is the possibility to exploit self-organizing phenomena for 
the fabrication of nanostructured surfaces by means of Molecular Beam 
Epitaxy (MBE) or similar techniques. 

In many cases, adsorbate and substrate materials crystallize in the same 
lattice but with different bulk lattice constants. Frequently, one observes 
that the adsorbate grows Jayer by layer, initially, with the lateral spacing of 
atoms adapted to the substrate. The misfit induces compressive or tensile 
strain in this pseudomorphic film and, eventually, misfit dislocations will 
appear. These relax the strain and the adsorbate grows with its natural 
lattice constant far from the substrate, eventually. 
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Dislocations will clearly dominate strain relaxation in sufficiently thick 
films and for large misfits. In material combinations with relatively small 
misfit an alternative effect governs the initial growth of thin films: Instead 
of growing layer by layer, the adsorbate aggregates in three-dim. structures. 
The term 3D-islands is commonly used to indicate that these structures are 
spatially separated. The situation is clearly different from the emergence of 
mounds due to the Ehrlich-Schwoebel (ES) instability [UlEI, for instance. 

At least two different growth scenarios display 3D-island formation: In 
Volmer-Weber growth, such structures appear immediately upon the sub- 
strate when depositing the mismatched material. The situation resembles 
the formation of non- wetting droplets of liquid on a surface. It is frequently 
observed in systems where adsorbate and substrate are fundamentally dif- 
ferent, an example being Pb on a graphite substrate HI. 

In the following we concentrate on the so-called Stranski-Krastanov 
(SK) growth mode, where 3D-islands are found upon a pseudomorphic 
wetting-layer (WL) of adsorbate material dElElE]. Two prominent pro- 
totype SK-systems are Ge/Si and InAs/GaAs where, as in almost all cases 
discussed in the literature, the adsorbate is under compression in the WL. 

In order to avoid conflicts with more detailed definitions and interpre- 
tations of the SK growth mode in the literature we will resort to the term 
SK-like growth in the following. It summarizes the following sequence of 
phenomena during the deposition of a few monolayers (ML) of material: 

1. The layer by layer growth of a pseudomorphic adsorbate WL up to 
a kinetic thickness h^ L . 

2. The sudden appearance of 3D-islands, marking the so-called 
2D-3D- or SK-transition 

3. Further growth of the 3D-islands, fed by additional deposition and 
by incorporation of surrounding WL atoms 

4. The observation of separated 3D-islands of similar shapes and sizes, 
on top of a WL with reduced stationary thickness hy/L- 

Besides these basic processes a variety of phenomena can play important 
roles in the SK-scenario, including the interdiffusion of materials and the 
segregation of compound adsorbates. These effects are certainly highly 
relevant in many cases, see f5j |6) and other contributions to this volume 
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iPfl . However, SK-like growth is observed in a variety of material systems 
which may or may not display these specific features. For instance, inter- 
mixing or segregation should be irrelevant in the somewhat exotic case of 
large organic molecules like PTCDA deposited on a metal substrate, e.g. 
Ag(lll). Nevertheless, this system displays SK-like growth in excellent 
accordance with the above operative definition 0. 

Despite the extensive investigation of SK-growth, a complete detailed 
theoretical picture is still lacking, apparently. This concerns in particular 
the nature of the 2D-3D transition. One problem clearly lies in the richness 
of the phenomenon. On the other hand, the very diversity of SK-like sys- 
tems gives rise to the hope that this growth scenario might be governed by 
a few basic universal mechanisms. Accordingly, it should be possible to 
capture and identify these essential features in relatively simple prototype 
systems. 

This hope motivates the investigation of simplifying models without 
aiming at the reproduction of material specific details. Some of the key 
questions in this context are: under which conditions does a WL emerge 
and persist? How does its thickness before and after the SK-transition 
depend on the growth conditions? Which microscopic processes trigger 
and control the sudden formation of 3D-islands? How do the island size 
and their spatial arrangement depend on the parameters of the system? 

Following earlier investigations of related phenomena, e.g. J9J [lOl HU 
[T2ll . we choose a classical pair potential ansatz to represent the interac- 
tions between atoms in our model. Here, we restrict our studies to the 
fairly simple case of a modified Lennard- Jones (LJ) system in 1 + 1 spatial 
dimensions, i.e. growth on a one-dim. substrate surface. As we interprete 
our model as a cross-section of the physical (2 + l)-dim. case, we still use 
the common term 2D-3D transition for the formation of multilayer from 
monolayer islands. 

We investigate our model by means of Kinetic Monte Carlo (KMC) 
simulation, see e.g. l[14ll for an introduction and overview. This concept 
has proven useful in the study of non-equilibrium dynamical systems in 
general and in particular in the context of epitaxial growth, see e.g. dE). 

Most frequently, pre-defined lattices are used for the representation of 
the crystal. So-called Solid-On-Solid (SOS) models which neglect lattice 
defects, bulk vacancies, or dislocations have been very successful in the in- 
vestigation of various relevant phenomena, including scenarios of kinetic 
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roughening or mound formation due to instabilities |T1|2|. There is, how- 
ever, no obvious way of including mismatch and strain effects in a lattice 
gas model. A potential route is to introduce additional elastic interactions 
between neighboring atoms in an effective fashion. In fact, such models 
of hetero-systems have been studied in some of the earliest Kinetic Monte 
Carlo simulations of epitaxial growth lfT31 fT6l . see ifTTl for a recent ex- 
ample of a so-called ball and spring model. In alternative approaches the 
strain field of a given configuration is evaluated using elasticity theory as, 
for instance, in lfl~8l . 

In order to account for strain effects more faithfully, including poten- 
tial deviations from a perfect lattice structure, it is essential to allow for 
continuous particle positions. Given at least an approximation for the in- 
teratomic potentials, a Molecular Dynamics (MD) type of simulation |[T9t 
would be clearly most realistic and desirable, see |20| for one example in 
the context of heteroepitaxial growth. However, this method suffers gen- 
erally from the restriction to short physical times on the order of 10~ 6 s or 
less. MBE relevant time scales of seconds or minutes do not seem feasible 
currently even when applying sophisticated acceleration techniques ||2T1 . 

Here, we put forward an off-lattice KMC method which has been in- 
troduced in JEHU [El and apply it in the context of the SK-scenario. Some 
of the results have been published previously in less detail [ 13 1. The pa- 
per is organized as follows: in the next section we outline the model and 
simulation method. Before analysing the actual SK-like growth in section 
|4]we present some basic results concerning various diffusion scenarios in 
section|3j In the last section we summarize and discuss open questions and 
potential extensions of our work. 



2 Model and method 

In our off-lattice model we consider pairwise interactions given by LJ- 
potentials of the form j\9\ 
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where the relative distance of particles % and j can vary continuously. 
As a widely used approximation we cut off interactions for > 3u. 
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The choice of the parameters {U a , cr} in Eq. ([TJ characterizes the dif- 
ferent material properties in our model: interactions between two substrate 
(adsorbate) particles are specified by the sets {U s ,a s = 1} and {U a ,a a }, 
respectively. Instead of a third independent set we set U as = y/U s U a , <J as = 
(a s + cr a ) /2 for the inter-species interaction. 

As the lattice constant of a monoatomic LJ crystal is proportional to cr, 
the relative misfit is given by e = (a a — a s ) /a s . Here we consider only 
cases with a a > a s , i.e. positive misfits e on the order of a few percent. 
For the simulation of SK-like growth we set U s > U as > U a . In such 
systems, the formation of a WL and potential layer by layer growth should 
be favorable, in principle. If not otherwise specified, we have set the misfit 
to e = 4%, a typical value for SK-systems, and used the LJ-prefactors 
U 8 = l.OeV, U a = 0.74eV, U as « 0.86eV. 

Growth takes place on a substrate represented by six atomic layers, 
with the bottom layer fixed and periodic horizontal boundary conditions. 
In the following we mostly refer to systems with L = 800 particles per 
substrate layer, additional simulations with L = 400 or 600 revealed no 
significant L-dependence of the results presented here. 

The deposition of single adsorbate particles is performed with a rate 
Rd = L F, where F is the deposition flux. As we interprete the substrate 
lattice constant as our unit of length, flux and deposition rate (measured in 
ML/s) assume the same numerical values. 

The rates of all other significant changes of the configuration are given 
by Arrhenius laws. We consider only hopping diffusion events at the sur- 
face and neglect bulk diffusion, exchange processes or other concerted 
moves. Furthermore, diffusion is restricted to adsorbate particles at the 
surface whereas jumps of substrate particles onto the surface are not con- 
sidered. As we will demonstrate in a forthcoming publication, these sim- 
plifications are justified for small misfits e in the LJ-system because the 
corresponding rates are extremely low, see also (HUD- 

The rate R4 for a particular event i is taken to be of the form 

Ri = v exp [—fif] (2) 

where T is the simulation temperature and k B the Boltzmann constant. For 
simplicity, we assume that the attempt frequency v is the same for all dif- 
fusion events. In order to relate to physical units we use v a = lO 12 ^ 1 
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wherever numerical results are given. The activation barriers E t are cal- 
culated on-line given the actual configuration of the system. This can be 
done by a minimal energy path saddle point calculation ifTOl \19\ . Here, 
we use a frozen crystal approximation which speeds up the calculation of 
barriers significantly, see J22J for an application in the context of strained 
surfaces. Note that the calculations are particularly simple in 1+1 dimen- 
sions: the path between neighboring local minima of the potential energy is 
uniquely determined and the transition state corresponds to the separating 
local maximum. 

An important modification concerns interlayer diffusion. LJ-systems 
in 1+1 dimensions display a strong additional barrier which hinders such 
moves at terrace edges (HE). This so-called Ehrlich-Schwoebel (ES) ef- 
fect is by far less pronounced in (2+l)-dim. systems, because interlayer 
moves follow a path through an energy saddle point rather than the pro- 
nounced maximum at the island edge. In our investigation of the SK-like 
scenario we remove the ES-barrier for all interlayer diffusion events by 
hand. One motivation is the above mentioned over-estimation in one di- 
mension. More importantly, we wish to investigate strain induced island 
formation without interference of the ES instability. Note that the latter 
leads to the formation of mounds even in homoepitaxy JUI3- 

The rates for deposition and diffusion are used in a rejection-free KMC 
simulation [fT4l . Using a binary search tree technique, one of the possible 
events j is drawn with probability Rj / R, where R = (Rd + J2iRi) is the 
total rate of all potential changes. Time is advanced by a random interval 
r according to the Poisson distribution P(r) = Re~ Rr ifLH . 

In order to avoid the artificial accumulation of strain due to inaccu- 
racies of the method, the potential energy should be taken to the nearest 
local minimum by variation of all particle positions in the system. This 
relaxation process affects both, adsorbate and substrate atoms. In order to 
reduce the computational effort, we restrict the variation to particles within 
a radius 3a s around the location of the latest event, in general. The global 
minimization procedure is performed only after a distinct number of steps. 
It is important to note, that both procedures do not lead to a substantial 
rearrangement. Significant changes of the activation energies due to relax- 
ation signal the necessity to perform the global procedure more frequently. 
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3 Diffusion processes 



Before analysing the SK-like scenario, we compare the barriers for hop- 
ping diffusion in various settings on the surface. The investigation of sys- 
tems like Ge/Si(001) reveal a very complicated scenario due to anisotropics 
and the influence of surface reconstructions For Ge on Ge(lll) 

the barrier for hopping diffusion is higher on the surface of a compressed 
crystal, whereas diffusion is faster on relaxed Ge ll24l . 

LJ- or similar models do not reproduce this feature, in general. 
Schroeder and Wolf ll22ll consider (2+l)-dim. single species LJ-sy stems 
and evaluate the diffusion barrier for a single adatom on surfaces in various 
lattice types. Among other results they find that mechanical compression 
of the crystal lowers the barrier for surface diffusion. However, in the mis- 
matched two species system, it is more important to compare diffusion on 
(a) the substrate, (b) the WL and (c) the surface of partially relaxed islands. 

The strong adsorbate-substrate interaction (U as > U a ) favors the for- 
mation of a WL, but it also yields deep energy minima and a relatively 
high diffusion barrier for adsorbate particles on the substrate. This effect is 
much weaker for particles on a complete wetting monolayer of adsorbate 
and, hence, the corresponding diffusion barrier is significantly lower. The 
faster diffusion further stabilizes the WL, as deposited particles will reach 
and fill in gaps easily. In principle, the trend extends to the following lay- 
ers. However, due to the short range nature of the LJ-potential the influence 
of the substrate essentially vanishes on WL of three or more monolayers. 
With the example choice e = 4% and U s = l.OeV, U a = 0.74eV, U as « 
0.866^ we find an activation barrier of E® ~ 0.57eV for adsorbate diffu- 
sion on the substrate and E\ ~ E% « 0A7eV for diffusion on the first and 
second adsorbate layer, respectively. 

For the SK-scenario the diffusion on islands of finite extension is par- 
ticularly relevant. Figure [TJ shows the barriers for diffusion hops on islands 
of various heights located upon a wetting monolayer. We wish to point 
at two important features: (1) Diffusion on top of islands is, in general, 
slower than on the WL and the difference increases with the island height. 
In our model, this is an effect of the partial relaxation or over-relaxation in 
the island top layer. (2) Depending on the lateral island size and its height, 
there is a more or less pronounced diffusion bias towards the island cen- 
ter, reflecting the spatially inhomogeneous relaxation. A similar effect has 
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Figure 1 : Diffusion barriers as obtained in our model for a single adatom 
on a flat symmetric multilayer island with 24 base particles and height 
1, 3, 5 layers (bottom to top curves). Symbols represent the activation en- 
ergies for hops from the particle position x to the left neighbor site. The 
island is placed on top of a single WL, interaction parameters are given in 
the text of section[3] The leftmost barriers correspond to downward jumps 
at the island edge with suppression of the ES-effect. Horizontal lines mark 
the barrier for adatom diffusion on the WL (lower line) and on perfectly 
relaxed adsorbate material (upper line). 



been observed in (2 + 1) dim. LJ-systems OS- 
Note that (2) has to be distinguished from the diffusion bias imposed 
by the ES-effect, which would be present even in homoepitaxy and with 
particle positions restricted to a perfect undisturbed lattice. 

Clearly, (1) and (2) favor the formation of islands upon islands and 
hence play an important role in SK-like growth. They concern adatoms 
which are deposited directly onto the islands as well as particles that hop 
upward at edges, potentially. As we will argue in the following section, up- 
ward diffusion moves play the more important role in the 2D-3D transition 
of our model. 



4 SK-like growth scenario 

In our investigation we follow a scenario which is frequently studied in 
experiments [1251 l26ll : In each simulation run a total of 4 ML adsorbate 
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Figure 2: A section of a simulated crystal as obtained for R d = 7.0ML/s 
and T = 500/^. Islands are located on a stationary WL with h c ~l, the six 
bottom layers represent the substrate. The darker a particle is displayed, 
the larger is the average distance from its nearest neighbors. 



material is deposited at rates in the range 0.5ML/s < Rd < 9.0ML/s. 
After deposition is complete, a relaxation period with Rd = of about 10 7 
diffusion steps follows, corresponding to a physical time on the order of 
0.3s. 

As we have demonstrated in [ 1 2 1 , strain relaxation through dislocations 
is not expected for, say, e = 4% within the first few adsorbate layers. In- 
deed misfit dislocations were observed in none of the simulations presented 
here. Results have been obtained on average over at least 15 independent 
simulation runs for each data point. 

In our simulations we observe the complete scenario of SK-like growth 
as described in the introduction. Illustrating mpeg movies of the simula- 
tions are available upon request or directly at our web pages lETl . A sec- 
tion of a simulated crystal after deposition and island formation is shown 
in FigureEl Ultimately, the formation of islands is driven by the relaxation 
of strain. As shown in the figure, material within the 3D-island and at its 
surface can assume a lattice constant close to that of bulk adsorbate. On 
the contrary, particles in the WL are forced to adapt the substrate structure. 

During deposition, monolayer islands located on the WL undergo a 
rapid transition to bilayer islands at a well-defined thickness /i WL - For the 
systematic determination of /i WL we follow [1251 and fit the density p of 3D- 
islands as p = p Q (h — /i WL ) a > finding comparable values of the exponent 
a. Figure E] displays the results for two different substrate temperatures 
T and various deposition rates Rd- The increase of /i WL with decreasing 
T agrees qualitatively with several experimental findings, see Eol as one 
example. Heyn discusses the effect of adsorbate/substrate interdiffusion 
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Figure 3: Kinetic critical WL thickness h^ L vs. the deposition flux for 
two different substrate temperatures. Both curves correspond to Eq. © 
with h and 7 obtained from the data of T = 500K, only. 

on the T-dependence of /i WL in the InGaAs system Q. Reassuringly, our 
result is consistent with his findings for absent intermixing. 

A key observation is the significant increase of /i WL with increasing 
deposition flux. Qualitatively the same flux dependence is reported for 
InP/GaAs heteroepitaxial growth in J26). This behavior leads to the con- 
clusion that the emergence of islands upon islands, i.e. the SK-transition, is 
mainly due to particles performing upward hops onto existing monolayer 
islands. If, on the contrary, the formation of second or third layer nuclei 
by freshly deposited adatoms was the dominant process, one would expect 
more frequent nucleation and an earlier 2D-3D transition at higher growth 
rates. Note again that we have suppressed the ES-effect explicitly which 
would also lead to more frequent mound formation at higher deposition 
rates Q0. 

In our model, upward diffusion is the limiting effect and sets the char- 
acteristic time for the SK-transition. The frequency of these processes 
strongly depends on the temperature and it has to be compared with the 
deposition rate. A high incoming flux will fill layers before particles can 
perform upward hops, and hence it will delay the 2D-3D-transition, cf. 
13 These considerations, together with the arguments of (26|, suggest a 
functional dependence of h^ L on F and T of the form 

h WL = h (^)\ (3) 
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Figure 4: Average base size (b) of multilayer islands as a function of 
Rd at T = 500fT, together with standard error bars. The inset shows the 
result for Rd = 4.5ML/s and different misfit parameters, the solid line 
corresponds to (b) = 0.91/e. 



where R up is the Arrhenius rate for upward hops. In our simulations we 
find a typical value which corresponds to an activation barrier E up rs 
1.0 eV close to the transition. Using this value we have performed a non- 
linear fit according to Eq. © based on the data for T = 500K. Inserting 
the obtained parameters h Q w 3.0ML and 7 ~ 0.2 in Eq. © yields good 
agreement for T = 480-fT, as well. Note that the value of 7 is expected to 
vary with the material system in the range < 7 < 0.5 according to ll26ll . 

Our assumption that upward diffusion is more important than nucle- 
ation of deposited adatoms on top of islands is further supported by the ob- 
servation that multilayer islands tend to form on a WL of, say, two mono- 
layers thickness even without particle deposition. 

After the 2D-3D-transition, islands grow by incorporating newly de- 
posited material, but also by consumption of the surrounding WL. Note 
that very large islands are observed to split by means of upward diffusion 
events onto their top layer E71 . The migration of WL particles towards 
and onto the islands can also extend into the relaxation period after depo- 
sition ends. Eventually, a stationary WL thickness /iwl ~ 1 is observed 
in our example scenario with U as ~ 0.86eV. By increasing the strength 
of the adsorbate/substrate interaction we can achieve, e.g., /i W l w 2 for 
U as ~ 2.7 'eV, but it is difficult to stabilize a greater stationary thickness, 
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due to the short range nature of the LJ-potential. 

Finally, we discuss the properties of 3D-islands emerging in SK- 
growth. As an example, their base length b is measured as the number of 
particles in the island bottom layer. The results discussed in the following 
were obtained at the end of the relaxation period with R d = 0. Whereas 
mean values do not change significantly, fluctuations are observed to de- 
crease with time in this phase. 

We observe that, for fixed temperature, the island size decreases with 
increasing deposition flux, as observed in several experimental studies 
1 26]. However, (b) becomes constant and independent of T for large 
enough deposition rate, cf. Figure |4j A corresponding behavior is found 
for the island density and their lateral distance, hinting at a considerable de- 
gree of spatial ordering [fl~3l . The saturation behavior further demonstrates 
the importance of upward hops vs. aggregation of deposited particles on 
islands. The latter process would yield a continuous increase of the island 
density with R d . 

We find a narrower distribution of island sizes with increasing R d lfT3l . 
The island size distribution in the saturation regime will be studied in 
greater detail in forthcoming investigations. 

We conclude our discussion by noting that, in the saturation regime of 
high growth rates, the typical island size follows a simple power law: (b) oc 
1/e, cf. Figure |4j Very far from equilibrium, the only relevant length scale 
in the system is given by the relative periodicity 1/e of the adsorbate and 
substrate lattices. This characteristic length was already found to dominate 
in the formation of misfit dislocations for larger values of e lfT2l . 

5 Conclusion and Outlook 

Despite its conceptual simplicity and the small number of free parameters, 
our model reproduces various phenomena of heteroepitaxial growth. We 
believe that, with a proper choice of interaction parameters, our model 
should be capable of reproducing all three prototype growth modes: layer 
by layer growth (for very small misfits), Volmer- Weber (for U as < U a ) and, 
as demonstrated here, the Stranski-Krastanov mode. Our work provides a 
fairly detailed and plausible picture of the latter. The key features are: 
(a) The strong adsorbate/substrate interaction favors the WL formation and 
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results in a relatively slow diffusion of adatoms on the substrate. Diffusion 
on the WL is significantly faster and the corresponding barrier decreases 
with the WL thickness. 

(b) Strain relaxation leads to a pronounced bias towards the island center on 
top of finite mono- and multilayer islands located on the WL. In addition, 
diffusion is slower on top of the partially relaxed islands than on the WL, 
in our model. 

Whereas (a) favors the formation and persistence of the WL, (b) clearly 
de-stabilizes layer by layer growth. We find that the microscopic process 
which triggers the transition is upward diffusion of adatoms from the WL 
and at island edges. The corresponding barriers decrease with the WL 
thickness analogous to (a). As a result of the competing effects, the 2D- 
3D-transition occurs at a critical thickness which depends upon T and R d 
as suggested in Eq. ©. 

We hypothesize that strain effects induce spatial modulations with a 
characteristic length scale e _1 and thus control the island size far from 
equilibrium. We find a corresponding saturation regime for large enough 
deposition fluxes. The precise mechanism of the size selection will be 
the subject of a forthcoming project. A related open question concerns 
the crossover from kinetically controlled island sizes to the equilibrium 
behavior which should be achieved in the limit of very small Rd. Several 
arguments [181 suggest that the typical island size close to equilibrium 
should be of order e~ 2 . 

Further investigations will concern the effects of intermixing and seg- 
regation which have been excluded from our model, so far. To this end 
we will consider the co-deposition of both species and allow for exchange 
diffusion at the substrate/adsorbate interface. 

In order to test the potential universality of our results, we will intro- 
duce different types of interaction potentials in our model. Ultimately, we 
plan to extend our model to the relevant case of 2 + 1 dimensions and to 
more realistic empirical potentials for semiconductor materials, e.g. ll28l . 
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